home *** CD-ROM | disk | FTP | other *** search
/ Technotools / Technotools (Chestnut CD-ROM)(1993).ISO / lang_oth / linpklib / sgeco.for < prev    next >
Text File  |  1983-12-31  |  6KB  |  194 lines
  1.       SUBROUTINE SGECO(A,LDA,N,IPVT,RCOND,Z)
  2.       INTEGER LDA,N,IPVT(1)
  3.       REAL A(LDA,1),Z(1)
  4.       REAL RCOND
  5. C
  6. C     SGECO FACTORS A REAL MATRIX BY GAUSSIAN ELIMINATION
  7. C     AND ESTIMATES THE CONDITION OF THE MATRIX.
  8. C
  9. C     IF  RCOND  IS NOT NEEDED, SGEFA IS SLIGHTLY FASTER.
  10. C     TO SOLVE  A*X = B , FOLLOW SGECO BY SGESL.
  11. C     TO COMPUTE  INVERSE(A)*C , FOLLOW SGECO BY SGESL.
  12. C     TO COMPUTE  DETERMINANT(A) , FOLLOW SGECO BY SGEDI.
  13. C     TO COMPUTE  INVERSE(A) , FOLLOW SGECO BY SGEDI.
  14. C
  15. C     ON ENTRY
  16. C
  17. C        A       REAL(LDA, N)
  18. C                THE MATRIX TO BE FACTORED.
  19. C
  20. C        LDA     INTEGER
  21. C                THE LEADING DIMENSION OF THE ARRAY  A .
  22. C
  23. C        N       INTEGER
  24. C                THE ORDER OF THE MATRIX  A .
  25. C
  26. C     ON RETURN
  27. C
  28. C        A       AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS
  29. C                WHICH WERE USED TO OBTAIN IT.
  30. C                THE FACTORIZATION CAN BE WRITTEN  A = L*U  WHERE
  31. C                L  IS A PRODUCT OF PERMUTATION AND UNIT LOWER
  32. C                TRIANGULAR MATRICES AND  U  IS UPPER TRIANGULAR.
  33. C
  34. C        IPVT    INTEGER(N)
  35. C                AN INTEGER VECTOR OF PIVOT INDICES.
  36. C
  37. C        RCOND   REAL
  38. C                AN ESTIMATE OF THE RECIPROCAL CONDITION OF  A .
  39. C                FOR THE SYSTEM  A*X = B , RELATIVE PERTURBATIONS
  40. C                IN  A  AND  B  OF SIZE  EPSILON  MAY CAUSE
  41. C                RELATIVE PERTURBATIONS IN  X  OF SIZE  EPSILON/RCOND .
  42. C                IF  RCOND  IS SO SMALL THAT THE LOGICAL EXPRESSION
  43. C                           1.0 + RCOND .EQ. 1.0
  44. C                IS TRUE, THEN  A  MAY BE SINGULAR TO WORKING
  45. C                PRECISION.  IN PARTICULAR,  RCOND  IS ZERO  IF
  46. C                EXACT SINGULARITY IS DETECTED OR THE ESTIMATE
  47. C                UNDERFLOWS.
  48. C
  49. C        Z       REAL(N)
  50. C                A WORK VECTOR WHOSE CONTENTS ARE USUALLY UNIMPORTANT.
  51. C                IF  A  IS CLOSE TO A SINGULAR MATRIX, THEN  Z  IS
  52. C                AN APPROXIMATE NULL VECTOR IN THE SENSE THAT
  53. C                NORM(A*Z) = RCOND*NORM(A)*NORM(Z) .
  54. C
  55. C     LINPACK. THIS VERSION DATED 08/14/78 .
  56. C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
  57. C
  58. C     SUBROUTINES AND FUNCTIONS
  59. C
  60. C     LINPACK SGEFA
  61. C     BLAS SAXPY,SDOT,SSCAL,SASUM
  62. C     FORTRAN ABS,AMAX1,SIGN
  63. C
  64. C     INTERNAL VARIABLES
  65. C
  66.       REAL SDOT,EK,T,WK,WKM
  67.       REAL ANORM,S,SASUM,SM,YNORM
  68.       INTEGER INFO,J,K,KB,KP1,L
  69. C
  70. C
  71. C     COMPUTE 1-NORM OF A
  72. C
  73.       ANORM = 0.0E0
  74.       DO 10 J = 1, N
  75.          ANORM = AMAX1(ANORM,SASUM(N,A(1,J),1))
  76.    10 CONTINUE
  77. C
  78. C     FACTOR
  79. C
  80.       CALL SGEFA(A,LDA,N,IPVT,INFO)
  81. C
  82. C     RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) .
  83. C     ESTIMATE = NORM(Z)/NORM(Y) WHERE  A*Z = Y  AND  TRANS(A)*Y = E .
  84. C     TRANS(A)  IS THE TRANSPOSE OF A .  THE COMPONENTS OF  E  ARE
  85. C     CHOSEN TO CAUSE MAXIMUM LOCAL GROWTH IN THE ELEMENTS OF W  WHERE
  86. C     TRANS(U)*W = E .  THE VECTORS ARE FREQUENTLY RESCALED TO AVOID
  87. C     OVERFLOW.
  88. C
  89. C     SOLVE TRANS(U)*W = E
  90. C
  91.       EK = 1.0E0
  92.       DO 20 J = 1, N
  93.          Z(J) = 0.0E0
  94.    20 CONTINUE
  95.       DO 100 K = 1, N
  96.          IF (Z(K) .NE. 0.0E0) EK = SIGN(EK,-Z(K))
  97.          IF (ABS(EK-Z(K)) .LE. ABS(A(K,K))) GO TO 30
  98.             S = ABS(A(K,K))/ABS(EK-Z(K))
  99.             CALL SSCAL(N,S,Z,1)
  100.             EK = S*EK
  101.    30    CONTINUE
  102.          WK = EK - Z(K)
  103.          WKM = -EK - Z(K)
  104.          S = ABS(WK)
  105.          SM = ABS(WKM)
  106.          IF (A(K,K) .EQ. 0.0E0) GO TO 40
  107.             WK = WK/A(K,K)
  108.             WKM = WKM/A(K,K)
  109.          GO TO 50
  110.    40    CONTINUE
  111.             WK = 1.0E0
  112.             WKM = 1.0E0
  113.    50    CONTINUE
  114.          KP1 = K + 1
  115.          IF (KP1 .GT. N) GO TO 90
  116.             DO 60 J = KP1, N
  117.                SM = SM + ABS(Z(J)+WKM*A(K,J))
  118.                Z(J) = Z(J) + WK*A(K,J)
  119.                S = S + ABS(Z(J))
  120.    60       CONTINUE
  121.             IF (S .GE. SM) GO TO 80
  122.                T = WKM - WK
  123.                WK = WKM
  124.                DO 70 J = KP1, N
  125.                   Z(J) = Z(J) + T*A(K,J)
  126.    70          CONTINUE
  127.    80       CONTINUE
  128.    90    CONTINUE
  129.          Z(K) = WK
  130.   100 CONTINUE
  131.       S = 1.0E0/SASUM(N,Z,1)
  132.       CALL SSCAL(N,S,Z,1)
  133. C
  134. C     SOLVE TRANS(L)*Y = W
  135. C
  136.       DO 120 KB = 1, N
  137.          K = N + 1 - KB
  138.          IF (K .LT. N) Z(K) = Z(K) + SDOT(N-K,A(K+1,K),1,Z(K+1),1)
  139.          IF (ABS(Z(K)) .LE. 1.0E0) GO TO 110
  140.             S = 1.0E0/ABS(Z(K))
  141.             CALL SSCAL(N,S,Z,1)
  142.   110    CONTINUE
  143.          L = IPVT(K)
  144.          T = Z(L)
  145.          Z(L) = Z(K)
  146.          Z(K) = T
  147.   120 CONTINUE
  148.       S = 1.0E0/SASUM(N,Z,1)
  149.       CALL SSCAL(N,S,Z,1)
  150. C
  151.       YNORM = 1.0E0
  152. C
  153. C     SOLVE L*V = Y
  154. C
  155.       DO 140 K = 1, N
  156.          L = IPVT(K)
  157.          T = Z(L)
  158.          Z(L) = Z(K)
  159.          Z(K) = T
  160.          IF (K .LT. N) CALL SAXPY(N-K,T,A(K+1,K),1,Z(K+1),1)
  161.          IF (ABS(Z(K)) .LE. 1.0E0) GO TO 130
  162.             S = 1.0E0/ABS(Z(K))
  163.             CALL SSCAL(N,S,Z,1)
  164.             YNORM = S*YNORM
  165.   130    CONTINUE
  166.   140 CONTINUE
  167.       S = 1.0E0/SASUM(N,Z,1)
  168.       CALL SSCAL(N,S,Z,1)
  169.       YNORM = S*YNORM
  170. C
  171. C     SOLVE  U*Z = V
  172. C
  173.       DO 160 KB = 1, N
  174.          K = N + 1 - KB
  175.          IF (ABS(Z(K)) .LE. ABS(A(K,K))) GO TO 150
  176.             S = ABS(A(K,K))/ABS(Z(K))
  177.             CALL SSCAL(N,S,Z,1)
  178.             YNORM = S*YNORM
  179.   150    CONTINUE
  180.          IF (A(K,K) .NE. 0.0E0) Z(K) = Z(K)/A(K,K)
  181.          IF (A(K,K) .EQ. 0.0E0) Z(K) = 1.0E0
  182.          T = -Z(K)
  183.          CALL SAXPY(K-1,T,A(1,K),1,Z(1),1)
  184.   160 CONTINUE
  185. C     MAKE ZNORM = 1.0
  186.       S = 1.0E0/SASUM(N,Z,1)
  187.       CALL SSCAL(N,S,Z,1)
  188.       YNORM = S*YNORM
  189. C
  190.       IF (ANORM .NE. 0.0E0) RCOND = YNORM/ANORM
  191.       IF (ANORM .EQ. 0.0E0) RCOND = 0.0E0
  192.       RETURN
  193.       END
  194.